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This  paper  proposes  a  probabilistic  approach  for  economic/emission  management  of  Micro  Grids 
(MGs).  In  order  to  meet  the  electrical  and  thermal  loads  while  having  lower  emission  production  in 
a  more  economical  manner,  combining  heating  and  power  along  with  hydrogen  production  and 
utilization  strategies  are  employed.  A  PEMFCPP  (Proton  Exchange  Membrane  Fuel  cell  power  plant)  is 
considered  as  the  prime  mover  of  the  Combined  Heat  and  Power  (CHP)  system.  The  surplus  power  of 
PEMFCPP  is  managed  to  produce  hydrogen.  An  electrochemical  model  for  representation  and 
performance  evaluation  of  the  PEMFC  is  applied.  Using  this  model,  the  output  voltage  and  power  of 
the  PEMFC  are  calculated  as  a  function  of  current,  constructive  and  operational  parameters.  The 
proposed  probabilistic  optimization  method  includes  2m  +  1  point  estimate  method  to  cover  the 
uncertainties  and  a  modified  multi-objective  algorithm  based  on  the  Modified  Gravitational  Search 
Algorithm  (MGSA)  to  find  Pareto-optimal  front  of  the  operation  management  problem.  The  study 
considers  the  uncertainties  in  forecasted  values  of:  the  hourly  market  tariffs,  electrical  and  thermal 
load  demands,  available  output  power  of  the  Photovoltaic  (PV)  and  Wind  Turbines  (WT)  units,  fuel 
prices,  hydrogen  selling  price,  operation  temperature  of  the  FC,  and  pressure  of  the  reactant  gases 
of  FC. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Micro  grids  are  clusters  of  distributed  energy  resource  (DER) 
units  located  near  the  electrical  and  thermal  loads,  and  can  be 
organized  to  work  in  both  grid-connected  and  autonomous  modes 
[1  ].  The  MGs  offer  a  wide  verity  of  benefits  to  the  utility  owners  and 
customers  such  as  [2]:  dealing  with  environmental  concern  by 
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reduction  of  carbon  emissions  because  of  using  renewable  energy 
sources  and  more  efficient  use  of  fossil  fuels,  energy  efficiency 
improvement  due  to  less  line  losses  and  co-generation  options, 
increasing  the  power  quality  and  reliability  because  of  generating 
power  near  the  loads,  finally  reducing  investment  risks  and 
increasing  fuel  diversification.  The  DER  units  can  be  categorized 
into  two  concepts:  the  distributed  generation  (DG)  and  the 
distributed  storage  (DS)  units.  The  DG  units  include  several  tech¬ 
nologies  such  as  FCs,  micro  turbines  (MTs),  diesel  engines,  PV  and 
WTs.  Likewise;  the  DS  units  contain  batteries,  energy  capacitors,  fly 
wheels  and  controllable  loads  [3]. 

Recently,  FCs  attract  much  attention  due  to  their  high  efficiency, 
low  aggression  to  the  environment,  vibration  free,  co-generation 
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Nomenclature 


Symbol,  description 

t,  k  time  interval  and  iteration  index,  respectively 
n  total  number  of  optimization  variables 

NT  total  number  of  hours 

Ng,  Ns  total  number  of  generation  expect  FCPP  and  storage 
units,  respectively 

Nd  total  number  of  load  levels 

/(X)  expected  cost 

u\  status  of  unit  i  at  hour  t 

,  pT  active  power  output  of  the  generator  i  and  the  storage 
device  j  at  time  t,  respectively 

PGrid  active  power  bought/sold  from/to  the  utility  at  time  t 

bgp  Bsj  bid  of  the  DG  source  i  and  the  storage  device  j  at  hour  t, 

respectively 

BGrid  bid  of  utility  at  hour  t 

SGi ,  Ssj  start-up/shut-down  costs  for  the  i  DG  unit  and  the 
storage  device  j,  respectively 
T\lect  total  electrical  load  at  time  t 
Pqi  min,  plQi  max  minimum  and  maximum  active  power 
production  of  the  DG  unit  i  at  hour  t 
Psj  min’  Psj  max  minimum  and  maximum  active  power  production 
of  the  storage  j  at  hour  t,  respectively 
Pgrid  min’  Pgrid  max  minimum  and  maximum  active  power 

production  of  the  utility  at  hour  t,  respectively 
Res f  the  scheduled  spinning  reserve  at  time  t 
P ceil, max  maximum  generating  power  limit  of  FCPP 
B£e//  price  of  natural  gas  for  FCPP 
B[h  fuel  price  for  residential  thermal  loads 

Bpump  hydrogen  pumping  cost 

B ^  hydrogen  selling  price 

P£e//  electrical  power  produced  by  FCPP  in  interval  t 
P^st  stored  hydrogen  amount  at  interval  t 
Pfj  equivalent  electric  power  for  hydrogen  production 
pH,end  available  amount  of  hydrogen  at  the  end  of  the  day 
Ph usage  secondary  hydrogen  stream  amount  in  kW  at  interval  t 

P^st  stored  hydrogen  amount  at  interval  t 

a ,  (3  hot  and  cold  start  up  cost  of  fuel  cell,  respectively 

electricalri  electrical  efficiency  at  full  load 
thermal tj  thermal  efficiency  at  full  load 
Reiect  resistance  to  the  electrons  flow 

RpWt  equivalent  membrane  resistivity 

Z/  value  of  the  Ith  input  random  variable  of  the  2m  point 

estimate  method 

zi,po  poth  standard  location  of  z* 

(x)ipo  poth  weighting  factor  z/>po 

ps,  as  mean  and  the  standard  deviation  of  the  S,  respectively 


N,  Prob (z;,j)  number  of  observations  of  Zi  and  the  probability  of 
each  observation  zij ,  respectively 
m  number  of  the  input  random  variables  of  2m  +  1  point 

estimate  method 

ftermax  maximum  number  of  iterations 

r  random  number  with  uniform  distribution  between 

0  and  1. 

Iter  current  iteration. 

R^e  Euclidian  distance  between  two  particles  j  and  e. 
e  small  positive  constant 

Nswrm  total  number  of  the  bees  in  the  swarm 
M|,m£  gravitational  mass  related  to  particles  j  and  e. 

Vnewj'Voidj  new  anc*  vel°city  of  the  jth  particle,  respectively. 

^newj^ofdj  new  anc*  old  position  of  the  jth  particle,  respectively 

Gk  gravitational  constant  at  the  kth  iteration 

Fit(Xj)  fitness  value  of  the  jth  particle. 

t  fuel  cell  cooling  time  constant 

rfcell  fuel  cell  efficiency  at  interval  t 

7]st  hydrogen  storage  efficiency 

l\h  thermal  load  demand  at  interval  t 

rTE  thermal  to  electrical  energy  ratio 

t0ff  time  the  FCPP  has  been  off 

iFC  cell  operating  current  (A) 

A  the  ratio  of  the  number  of  water  moles  for  each 

sulfonic  group  in  the  membrane 
T  cell  operation  temperature  (I<) 

List  of  abbreviations 
MG  micro  grid 

MGSA  modified  gravitational  search  algorithm 

GSA  gravitational  search  algorithm 

DER  distributed  energy  resource 

DG  distributed  generation 

DS  distributed  storage 

WT  wind  turbine 

MT  micro  turbine 

PV  photovoltaic 

EEM  economic/emission  management 

MGCC  micro  grid  central  controller 

PEMFC  proton  exchange  membrane  fuel  cell 

FCPP  fuel  cell  power  plant 

NiMFI-battery  nickel-metal-hydride  battery 

PDF  probability  density  function 

CDF  cumulative  density  function 

IRV  input  random  variable 

ORV  output  random  variable 

DM  decision  maker 

STD  standard  deviation 


options,  superior  reliability  and  modular  nature  [4].  In  particular, 
PEMFCs  emerge  as  a  promising  candidate  for  both  stationary  and 
automotive  applications.  The  attractiveness  of  the  PEMFCs  in 
automotive  applications  is  due  to  their  high  power  density  (about 
0.5  W/cm2)  and  low  operating  temperature  (about  50-90  °C).  The 
PEMFCs  consume  pure  hydrogen  as  fuel  and  produce  water  as  a  by¬ 
product  waste  [5].  Typically  one  single  cell  produces  about  0.5- 
0.9  V,  so  in  power  generating  systems  several  cells  are  connected 
in  series  to  create  a  stack  which  can  provide  hundreds  of  kilowatt 
[6].  Several  studies  have  shown  that  in  order  to  make  the  FCs  more 


cost-effective  it  is  necessary  to  recover  waste  thermal  energy  and 
manage  producing  hydrogen  [7,8]. 

Flydrogen  is  one  of  the  most  environmental  friendly  and  bene¬ 
ficial  fuels  which  can  operate  as  a  storage  medium  or  energy  carrier. 
Reaction  of  hydrogen  with  oxygen  produces  water  with  no 
pollutant  [9-11  ].  Generally,  the  production  of  hydrogen  from  fossil 
fuels  is  considered  on  an  ‘as  needed’  basis.  Flowever,  sometimes  it  is 
more  economical  and  efficient  to  store  the  hydrogen  fuel  as 
hydrogen  [7].  In  this  regard,  in  order  to  make  FCPPs  more  profitable, 
some  amount  of  the  reformer  capacity  can  be  employed  to  produce 
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hydrogen  to  be  stored.  The  stored  hydrogen  can  be  reused  to 
generate  electricity  by  the  FCs  or  to  be  sold  to  other  customers. 
Moreover,  recovering  the  waste  thermal  energy  of  the  FCs  to  make 
a  CFIP  system  is  one  of  the  significant  options  to  improve  the  FCs 
efficiency.  More  efficient  energy  used  in  the  CHP  systems  leads  to 
cost  and  fuel  saving  and  consequently  carbon  emission  reduction 
[12]. 

In  the  last  few  years,  several  researches  have  focused  on 
environmental  and  economic  concepts  of  the  MGs,  the  CFIP 
systems  and  hydrogen  production.  Chen  et  al.  [3]  presented 
a  smart  energy  management  system  based  on  the  matrix  real- 
coded  genetic  algorithm  to  optimize  the  operation  of  the  MGs. 
The  proposed  framework  in  Ref.  [3]  consists  of  a  power  forecasting 
module,  an  energy  storage  system  management  module  and  an 
optimization  module.  Flernandez  et  al.  developed  an  optimization 
scheme  to  reduce  the  fuel  consumption  of  the  MGs  while  the  local 
energy  demand  (both  electrical  and  thermal)  and  a  certain 
minimum  reserve  power  were  satisfied.  Tsikalakis  et  al.  optimized 
the  operation  of  the  MGs  during  inter-connected  operation  by 
optimizing  the  generation  of  the  DG  units  and  power  exchange 
with  the  upstream  network  [13].  In  Ref.  [14],  a  linear  programming 
was  examined  to  minimize  the  average  production  cost  of  electric 
power  in  a  hybrid  solar-wind  MG  while  environmental  factors 
were  considered.  Moreover,  in  the  filed  of  the  FC  micro  CFIP  and 
hydrogen  management,  El-Sharkh  et  al.  proposed  an  economic 
model  for  PEMFC  and  used  an  evolutionary-based  technique  to 
find  an  operational  strategy  for  hydrogen  utilization  and  waste 
thermal  energy  recovery  [15].  The  introduced  model  of  the  PEMFC 
in  Ref.  [15]  is  an  imperfect  one  which  doesn’t  consider  the  oper¬ 
ation  condition  of  the  FC  such  as  its  operation  temperature  and 
pressure  of  reactant  gases.  Moreover,  the  aforementioned  paper 
doesn’t  analyze  environmental  criterion.  In  Ref.  [16],  an  experi¬ 
mental  investigation  of  the  PEMFC  for  both  heat  and  power 
production  was  discussed  and  shown  that  utilization  of  the 
generated  heat  for  domestic  water  heating  could  increase  the 
efficiency  of  the  FCs  to  about  70%  in  comparison  with  only  35-50% 
for  electricity  generation  alone.  In  Ref.  [17],  an  economic  and 
environmental  evaluation  of  two  micro  CFIP  alternatives,  namely, 
gas  engines  and  FC  with  different  operating  modes  for  residential 
buildings  were  proposed.  The  analyses  demonstrate  that  the  FC 
systems  could  select  as  a  better  options  for  the  examined  resi¬ 
dential  building  from  both  environmental  and  economic  points  of 
view.  Hawkes  et  al.  devised  a  techno-economic  model  for  the 
design  and  control  of  FCs  micro  CHP  system  and  discussed  the 
impact  of  the  stack  degradation  on  economic  and  environmental 
performances  [18,19].  In  Ref.  [20],  a  comparison  between  the  heat 
led  and  electricity  led  operation  strategies  for  a  residential  micro 
CHP  system  was  offered. 

One  of  the  drawbacks  associated  with  previous  studies  for  the 
operation  management  of  the  MGs  is  neglecting  the  uncertainty  in 
generation  patterns,  load  demand  and  market  tariffs,  fuel  prices 
and  operation  condition  of  the  FCs.  The  deterministic  approaches 
are  dependent  on  the  accuracy  of  the  input  data,  so  errors  in  the 
input  random  data  leads  to  unreliable  solutions.  In  deed,  some 
phenomena  have  uncontrollable  nature  or  too  implicit  circum¬ 
stance  for  being  accurately  modeled.  The  market  tariffs,  fuel  prices 
and  load  demands,  for  example,  regarding  to  open  access  market 
and  divers  customer  types  are  more  unpredictable  than  before  [21- 
23].  Moreover,  due  to  the  stochastic  nature  of  the  wind  speed  and 
the  solar  radiation,  the  PV  and  WT  units  generate  uncontrollable 
and  fluctuated  power.  Therefore,  the  validity  of  the  traditional 
optimization  methods  should  be  re-examined  under  new  circum¬ 
stance  [13].  In  this  regard,  from  the  system  planning  point  of  view, 
new  approaches  need  to  be  employed  for  the  economic/emission 


management  (EEM)  of  the  MGs  to  cope  with  the  intermittency  in 
input  random  data  and  minimize  the  risk  associated  with  the 
design  and  operation  under  uncertainty  [24,25].  This  paper  for 
handling  the  uncertainty  organizes  the  EEM  problem  as  a  probabi¬ 
listic  one. 

The  probabilistic  methods  can  be  classified  in  three  categories: 
the  Monte  Carlo  Simulation  (MCS)  [26],  the  analytical  techniques 
and  the  approximate  methods.  The  MCS  randomly  generates 
several  values  for  each  input  random  variable  (IRV)  and  solves  the 
problem  for  each  of  them  as  a  deterministic  one.  Although  MCS 
can  provide  accurate  solutions,  its  huge  computational  burden 
makes  it  unattractive  and  unacceptable  for  practical  applications. 
Analytical  techniques  apply  fewer  numbers  of  simulations  but 
require  more  assumptions  and  complicated  mathematical 
computations  [27].  Approximate  methods  provide  a  good  balance 
between  computational  efficiency  and  accuracy  [28].  The  Point 
estimate  method,  as  an  approximate  method,  is  an  efficient  and 
reliable  method  to  model  the  uncertainty  in  power  systems  [19]. 
Mostly,  the  2m  point  estimate  method  has  been  used  to  deal  with 
power  system  problems.  The  2m  scheme  generally  fails  to  give 
satisfactory  results  when  the  number  of  uncertain  input  variables 
is  high  [25]. 

In  this  paper,  the  2m  +  1  point  estimate  method  is  devised  to 
optimize  the  energy  management  of  MGs.  The  estimate  points,  in 
this  method,  are  independence  to  the  number  of  IRVs  thus  the 
2m  +  1  scheme  unlike  the  2m  one  can  provide  satisfactory  results 
in  problems  with  a  large  numbers  of  uncertain  inputs  [25].  The 
2m  +  1  scheme  generates  three  values  for  each  IRV  and  like  the 
MSC  solves  the  problem  by  deterministic  routines  for  them.  In  this 
study,  the  hourly  market  tariffs,  electrical  and  thermal  load 
demand,  available  output  power  of  the  PV  and  WT  units,  fuel 
prices,  operation  temperature  of  the  FC,  pressure  of  the  reactant 
gases  of  FC  are  considered  as  IRVs.  The  normal  and  beta  proba¬ 
bility  density  function  (PDF)  are  used  to  model  the  variations  of 
IRVs.  Moreover,  the  Gram-Charlier  expansion  is  employed  to 
provide  more  accurate  probability  distribution  for  Output  Random 
Variables  (ORVs).  An  electrochemical-based  FC  model  is  used  to 
calculate  output  power  of  the  FC  as  a  function  of  load  current  and 
membrane  humidity.  Furthermore,  total  operational  cost  and 
emissions  of  the  MGs  including  emission  and  power  production  of 
different  DG  and  DS  units,  trade  with  the  market  and  hydrogen 
utilization  management  are  formulated.  The  EEM  model  empha¬ 
sizes  a  practical  formulation  of  the  PEMFC  and  links  it  to  the 
economic  and  emission  model.  Finally,  a  multi-objective  modified 
gravitational  search  algorithm  (MGSA)  is  devised  to  optimize  the 
EEM  of  the  MGs. 

The  GSA  is  inspired  by  the  law  of  gravity  and  mass  interac¬ 
tions.  In  GSA  method,  the  individuals  are  collection  of  masses 
which  interact  with  each  other  using  a  policy  inspired  by  the 
Newtonian  gravity  law  and  also  the  laws  of  motion  [29].  This 
study  takes  advantages  of  the  GSA  and  improves  it  to  devise 
a  powerful  multi-objective  optimization  algorithm  for  solving  the 
probabilistic  EEM  problem  of  the  MGs.  Firstly,  in  order  to  enhance 
the  diversity  of  solutions,  a  mutation  technique  is  suggested.  It  is 
expected  this  mutation  to  help  the  algorithm  to  escape  being 
trapped  into  local  optima.  After  that,  the  modified  GSA  is 
expanded  to  make  a  multi-objective  one.  Finally,  the  multi¬ 
objective  MGSA  is  linked  to  2m  +  1  point  estimate  method  to 
find  the  uniformly  POF  of  the  probabilistic  EEM  problem  with  two 
conflicting  objectives  (cost  and  emission  minimization).  The 
effectiveness  of  the  proposed  technique  is  verified  on  a  typical 
MG  participates  in  the  open  market  considering  the  FC,  electrical 
and  hydrogen  storage,  MT,  PV  and  WT  units  and  hydrogen 
production. 
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2.  Operation  management  of  a  micro  grid 

2.1.  Objective  functions 

2  A  A.  Minimization  of  the  total  operating  cost 

The  total  operating  cost  of  the  MG  includes  the  fuel  costs  of  units 
as  well  as  their  start-up/shut-down  costs.  The  objective  function 
can  be  formulated  as  follows  [15,30]: 

NT  NT  (  Ng 

Minf 1  (X)  =  J2Costt  =  J2 \uiPaBa  +Startci 

t=  1  t=  1  {  2  =  1 

x  max^O ,u- +5/iutQ-  x  max(o,u-_1  -ufj  j 

Ns 

+  ^  + Start Sj  X  max (0,  uj  -  uj”1  J  +  Shutsj 

j= i 

X  max(o, uj-1  -  uj)]  +PcCridBcCrid 
+  5ce/Z  (  (Pcell  +  Ph)  / Pcell)  Bth  x  max  (Lth  ~  Pth 1  °) 

BpumpVstBH  |  +  ^  +  /2  ^1  —  exp  (toff  A)  )  —  BHsPn.end 

+  OM  (1) 

where  X  =  [X^.-.X^.X^]  and  Xf  is  design  variables  vector 
including  active  powers  of  units  and  their  related  states  that  can  be 
described  as  follows: 


2.2.  Constraints 

•  Power  balance: 

Ng  Ns 

E^Gi  +  £p„  +  PGrid  +  Pcell  +  PH -usage  =  ^ elect 
i=l  j=l 

•  Real  power  generation  capacity: 

^Gi, min  —  ^Gi  —  Pd, max 
Psi,  min  —  Psi  —  Psi,  max 

(6) 

Dt  <  Pt  <  Pf 

grid,  min  —  Grid  —  grid,  max 

^ce/Z, min  <4«  +  pu  —  Pcell,m  ax 

•  Spinning  reserve: 

Ng  Ns 

X>iPci, max +  Z  Psj, max  +  P grid, max  + Pcell, max  ^  LeZect  +  ( 7 ) 

i=l  j=l 

3.  CHP  system 

In  this  study,  PEMFC  is  considered  as  a  prime  mover  of  CHP 
system.  The  PEMFC  consumes  pure  hydrogen  as  its  fuel;  transform 


Xf  = 


[pci  5 Pg2 j  •  •  -5 PGNg >  ^Grid ’  ^ceZZ ’PH’Pth’Psl’Ps2’  •**» PsNs >  U1  j  u2 >  ■  •  •  j  uNs+Ng+2] 


(2) 


2.1.2.  Minimization  of  the  total  pollutants  emissions 

In  this  study,  the  environmental  effects  of  carbon  dioxide  (CO2), 
sulfur  dioxide  (SO2)  and  nitrogen  oxides  (NOx)  are  involved  in  the 
optimization  problem.  Mathematical  functions  which  formulate 
emissions  associated  with  the  power  production  of  different  DG 
units  are  employed.  The  second  objective  can  be  described  as 
follows: 


Ns 


Minf2(X)  =  J2  Emission1  =  jT  A  Pc,  A]  +  £  WjP'sAi 

t=l  t=l  [  1  =  1  j= 1 

+  PGridEGrid  \  (3) 


where  E^,  and  *Crid  define  the  amount  of  emissions  in  kg/kWh 
for  each  DG,  storage  unit  and  the  utility  at  hour  t  respectively.  These 
variables  are  expressed  as  follows: 


Ef  • 

unite 


=  COrn 


-SOVn 


NOt, 


(4) 


where  CO?  ,  SO?  ,NOl  and  E\-,  are  the  amounts  of  carbon 

unite  unite  ^unite  UlllLe 

dioxide,  sulfur  dioxide  and  nitrogen  oxides  emissions  from  unit  e 
(DG/DS/utility)  at  hour  t,  respectively. 


the  chemical  energy  liberated  during  the  electrochemical  reaction 
of  hydrogen  and  oxygen  to  electrical  energy  and  produce  zero 
emissions  besides  water  vapor.  Typical  configuration  of  a  PEMFC 
can  be  found  in  Fig.  1  [5].  Separation  of  hydrogen  from 
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Fig.  1.  Basic  PEMFC  operation. 
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hydrocarbons  is  carried  out  through  an  exothermic  reaction  known 
as  reforming  and  is  done  by  reformers.  It’s  assumed  that  some 
amount  of  produced  hydrogen  by  reformer  is  stored  for  further 
usage  also  a  part  of  local  thermal  loads  are  satisfied  by  recovered 
thermal  energy  of  the  reformer. 

3.1.  Hydrogen  production  strategy 

As  shown  in  Fig.  2,  the  generated  power  of  the  reformer  divided 
into  two  parts.  One  of  them,  by  the  main  hydrogen  stream  enters  to 
the  stack;  while  the  other  part  is  stored  in  the  hydrogen  tank  for 
further  utilization.  In  order  to  take  account  of  the  later  part  of  the 
generated  hydrogen  in  the  PEMFC  model  an  equivalent  electrical 
power  is  attributed  to  it  as PH  [15].  In  this  study,  Ph  is  considered  as 
one  of  the  design  variables  and  its  value  can  change  in  the  range  of 
zero  and  Ppiff,  where  Pm/ji s  defined  as: 

P Diff  =  Pcell,  max  —  ^cell  (8) 

As  in  this  paper,  the  recovered  energy  of  the  reformer  is  used  for 
co-generating,  in  high  thermal  demand  intervals  which  reformer 
works  in  high  capacity  to  satisfy  thermal  load  demands  as  much  as 
is  economical,  some  amount  of  the  generated  hydrogen  are  stored 
in  the  hydrogen  tank.  The  stored  hydrogen  in  addition  to  the 
hydrogen  produced  by  the  reformer  can  be  used  in  the  low  thermal 
demand  intervals  to  generate  the  required  electrical  power.  The 
unused  hydrogen  in  the  hydrogen  tank  is  sold  at  the  end  of  the  day. 
The  hydrogen  production  in  kg/s  can  be  computed  as  follows  [15]: 

H 2(Amount)  =  1.05  x  108  x  (yCj  (9) 

The  stored  hydrogen  amount  in  the  hydrogen  tank  can  be 
calculated  as  follows: 

PHst  =  PHst  +PHX  Vst  ~  PHusage  (1°) 

It  is  anticipated  that  this  policy  makes  the  operation  manage¬ 
ment  of  the  MG  more  economical  and  improves  the  efficiency  of 


the  FCPP.  In  this  paper,  the  hydrogen  storage  cost  is  assumed  equal 
to  the  pumping  cost  and  no  storage  infrastructure  or  technology 
cost  is  taken  to  account. 

3.2.  Recovered  thermal  energy  FCPP 

This  work  recovers  the  waste  heat  generated  by  the  reformer 
in  order  to  make  a  CFIP  system  and  satisfy  some  amount  of  the 
MG  thermal  loads  by  recovered  thermal  energy.  It  is  expected 
that  the  CHP  systems  provide  a  higher  thermal  efficiency  in 
comparison  with  producing  heat  and  power  separately;  conse¬ 
quently  for  the  same  output,  less  fuel  is  consumed  which  causes 
to  reduce  the  operational  cost  and  emissions  of  greenhouse  gases 
[31].  The  recovered  thermal  power  from  the  FCPP  can  be 
computed  as: 

P[h  =  rrex(i^  +  i^B)  (11) 

where  rTE  along  with  the  efficiency  of  the  cell  rfcell  are  function  of 
PLR  (PLR  =  electrical  generated  power/maximum  power)  [15]. 
Recovering  the  thermal  power  and  the  hydrogen  management 
improve  the  overall  efficiency  of  the  cell  greatly.  The  overall  effi¬ 
ciency  can  be  given  as  [15]: 

„t  =  PceH  +  ??stxPH  +  m 
I overall  pt  ,  pt 

1  cell  +  1  H 

rfcell 

Moreover,  Thermal  and  electrical  efficiencies  at  full  load  can  be 
defined  by  the  following  formula  [32,33]: 

electrical ij  =  Pcel1  st  x  Ph  (13) 

'  in 

thermal p  =  m^n^hihhl  (14) 

l  in 


[  J  ]  Pressure  regular 


DI  water 


Water  pump 


Fig.  2.  Generation  system  with  PEMFC  stacks. 
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4.  2m  +  1  point  estimate  method 

This  paper  employs  2m  +  1  point  estimate  method  to  cover  the 
uncertainty  in  forecasted  values  of  electrical  and  thermal  load 
demands,  available  output  power  of  the  PV  and  WT  units,  market 
prices  in  each  interval,  natural  gas  price  for  the  FCPP,  operation 
temperature  of  the  FC,  pressure  of  the  reactant  gases  of  FC,  fuel 
price  for  those  thermal  loads  of  the  MG  which  are  not  satisfied  by 
the  FCPP  and  hydrogen  selling  price. 

Firstly,  the  point  estimate  method  was  introduced  by  Rose- 
nblueth  in  1975  [34]  but  there  is  a  significant  drawback  associated 
with  Rosenblueth’s  method  and  that  is  their  heavy  computational 
burden  requirement.  In  1989,  Harr  developed  a  new  point  estimate 
method  to  overcome  the  drawback  of  Rosenblueth’s  method  [35]. 
Although  the  Harr’s  method  was  computationally  more  efficient 
than  the  prior  ones,  it  was  restricted  to  symmetric  variables.  In 
1998,  Hong  introduced  a  new  point  estimate  method  which  was 
computational  more  efficient  than  the  previous  ones  [36].  The 
Hong’s  method  is  suitable  for  both  symmetric  and  asymmetric 
variables.  This  paper  proposes  the  Hong’s  2m  +  1  point  estimate 
method  to  solve  the  probabilistic  EEM  problem  of  the  MGs.  Mostly, 
the  Hong’s  2m  point  estimate  method  has  been  used  in  the  litera¬ 
ture  for  power  system  problems.  The  2m  +  1  scheme  is  more 
accurate  than  the  2m  one  for  two  reasons.  First,  the  2m  +  1  scheme 
considers  both  skewness  and  kurtosis  of  IRVs  while  the  2m  scheme 
doesn’t  use  the  second  one.  In  addition,  the  estimated  points  in  the 
2m  +  1  scheme  unlike  those  in  the  2m  one  are  independent  to  the 
number  of  IRVs,  so  in  problems  with  a  large  numbers  of  uncertain 
variables  again  it  can  be  an  efficient  method.  It  is  worth  mentioning 
that,  the  2m  +  1  scheme  needs  only  one  additional  evaluation  of  the 
objective  function  in  composition  with  the  2m  scheme. 

The  Deterministic  optimal  planning  of  the  MGs  depends  on  a  set 
variables  v  such  as  output  power  of  different  units,  load  demands 
and  market  prices  as: 

Sorv  =f(v)  (15) 

where  Sorv(ORV  =  1,2)  are  the  outputs  of  optimal  planning  and  / 
is  the  set  of  operation  cost  and  emission  equations. 

In  the  deterministic  EEM  of  an  MG,  the  aforementioned  vari¬ 
ables  have  fixed  values  but  in  realty,  some  of  those  variables  are 
uncertain.  For  example,  there  are  always  errors  in  forecasted  values 
for  the  output  power  of  the  WTs,  the  PVs  or  the  market  price  and 
the  load  demands.  The  function  /transfers  the  uncertainty  from  the 
IRVs  to  the  ORVs  as  follows: 

$orv  =  /(G ,  z2 , . . . ,  zm)  (16) 

where  c  is  the  set  of  certain  variables,  Zi(l  =  l,...,m)  are  input 
variables  under  uncertainty  and  Sorv  are  the  outputs  of  optimal 
planning. 

The  goal  of  the  2m  +  1  scheme  is  finding  the  statistical  infor¬ 
mation  of  the  ORVs  by  using  the  first  few  central  moments  of  the 
IRVs  i.e.  the  mean  pP(,  variance  oPl,  skewness  XPl  3  and  kurtosis  XPl4 
coefficients.  Employing  just  first  few  central  moments  of  IRVs  to 
evaluate  the  characteristics  of  the  output  set  is  a  remarkable 
advantage  of  the  point  estimate  methods  where  implementing  the 
features  of  IRVs  is  difficult  task  to  reach. 

The2m  +  1  method  produces  two  probability  concentration  for 
each  IRV  zi  as  (z/,i,Wji)  and(z/2,  w/)2).  Thez/po(po  =  1,2)  is  called 
the  poth  location  of  z/  and  w/po(po  =  1,2)  is  a  weighting  factor 
which  specifies  the  importance  of  the  corresponding  location  in 
evaluating  the  statistical  moments  of  the  ORVs.  The  deterministic 
EEM  problem  is  simulated  2m  +  1  times  in  the  proposed  probabi¬ 
listic  method.  One  of  the  mentioned  simulations  is  done  by  fixing 


all  the  IRVs  on  their  mean  values  called  Mean  location  while  in  each 
of  the  2m  reaming  simulations,  one  of  the  IRV  is  fixed  to  one  of  its 
locations,  and  the  other  IRVs  are  equal  to  their  mean  values  as 
follows: 

SoRV(i,po)  ~/^GMzi  ’^Z2’  5  po  =  1,2  l  =  1,2,  ...,m 

(17) 

where  c  is  the  set  of  certain  variables,  z/ti  and  zi  2  are  the  specified 
locations  of  the  IRV  zj  and  pZ)  is  the  mean  value  of  the  left  over  IRVs. 
Once  the  solutions  of  2m  +  1  deterministic  EEM  are  explored,  using 
the  weighting  factors  the  first  and  the  second  moments  of  each 
ORVs  can  be  estimated. 

The  step  by  step  procedure  of  2m  +  1  point  estimate  method  to 
calculate  the  moments  of  the  ORVs  can  be  summarized  as 

Step  1 :  Definem. 

Step  2:  Set  E(S£W)  =  0,  h  =  1,  2. 

Step  3:  Select  an  uncertain  parameterz/. 

Step  4:  Calculate  the  skewness  (AZ)3)  and  the  kurtosis  (XZl4)  of 
the  z/  according  to  the  following  equations: 


(19) 


Step  5:  Calculate  two  standard  locations: 


?tpo  =  ^f~  +  (-l)3-p°yAz„4-|43,  po  =  1,2  (20) 

Step  6:  Compute  two  estimated  location: 


Zf.po  =  Mz,  +  ?z,,po-<7z„  po  =  1,2  (21) 

Step  7:  Calculate  the  deterministic  EEM  for  the  poth  estimated 
location: 


SoRV(i,po)  —  ?Mz2 5 •  •  po  —  1,2  l  —  l,2,...,m 

(22) 

Step  8:  Compute  two  weighting  factors  of  zp 


w,,po  =  (-l)3-p7?/,po(?u  -fw).  VO  =  1,2  (23) 

Step  9:  Update  the  first  and  second  moment  of  the  ORVs  (total 
operational  cost  and  emission): 
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:($orv)  —  ^(Sorv)  +  w/,po '  ($ORV(l,po))  ,  h  —  1,2  (24) 

po  =  1 


Step  10:  Repeat  steps  3-9  until  all  uncertain  parameters  were 
taken  into  account. 

Step  11:  Compute  the  weight  factor  of  the  mean  location  as 
follows: 


(25) 

Step  12:  Calculate  the  deterministic  EEM  for  the  mean  location 
as: 


Sorv,h  =  l  =  1,2,  ...,m  (26) 


G*,G„xexp(s,>!jj|?-)  (30) 

where  G0  and  rn  are  two  constant  which  are  set  to  100  and  20, 
respectively  [29]. 

The  value  of  mass  of  each  particle  is  computed  by  mapping  its 
fitness  value  as  the  following  equations: 

fit(xf)-  worst* 

m‘~  best-  worst*  (3,) 


M'l 


rrij 

ST^Nswrm 

2^e  =  l 


(32) 


where  worstk  and  best k  are  the  maximum  and  the  minimum  fitness 
values  at  iteration  k  (in  minimization  problems). 

The  overall  gravitational  force  exerted  on  the  jth  particle  is 
calculated  by  a  randomly  weighted  sum  of  forces  applied  by  the 
other  particles  as  follows: 


Step  13:  compute  the  first  and  second  moment  of  the  ORV  using 
(Ofi  andSM: 


E(sh0RV)  =  E(shORV)+M^(SoRV^h,  h  =  1,2  (27) 


=  E  rex  Iff 

e  =  1  ,e^j 


(33) 


The  random  parameter  re  is  added  to  (33)  to  give  the  algorithm 
stochastic  characteristics.  The  acceleration  of  each  particle  is 
defined  as  follows: 


Step  14:  Compute  the  mean  and  standard  deviation  of  the  total 
operational  cost  and  emission. 


( SORV )  (^(5(W)) 


2 


(28) 


•k.t 


F 

aks  =  J- 

J  Mi 


(34) 


Finally,  the  new  position  and  velocity  of  each  particle  are 
updated  by  the  following  equations: 


yk,t  _  k,t  .  yk,t 
vnewj  uj  +  v old  j 


(35) 


The  probability  density  function  of  each  ORV  can  be  approxi¬ 
mated  and  plotted  using  its  calculated  mean  and  standard  devia¬ 
tion  and  also  Gram-Charlier  series  approach  [37]. 

5.  Modified  GSA 

5.1.  Overview  of  standard  GSA 


w7<,t  _  w7<,t  -I- 

Anew.j  —  Ao/dj  '  V new  j 


(36) 


In  order  to  balance  the  exploration  and  exploitation  capability  of 
the  GSA,  only  a  set  of  the  particles  with  better  fitness  values,  i.e. 
bigger  mass,  are  selected  to  apply  their  forces  to  the  others.  Hence, 
(33)  can  be  rewritten  as: 


The  GSA  has  inspired  by  the  Newton’s  law  of  gravity  introduced 
in  2009  [29].  The  gravitation,  as  one  of  the  fundamental  interac¬ 
tions  in  the  nature,  is  the  tendency  of  objects  to  accelerate  toward 
each  other.  According  to  the  law  of  gravity,  each  mass  attracts  every 
other  one  with  a  ‘gravitational  force’.  The  gravitational  force 
between  two  particles  is  directly  proportional  to  the  product  of 
their  masses  and  inversely  proportional  to  the  square  of  the 
distance  between  them.  Moreover,  based  on  the  Newton’s  second 
law,  the  acceleration  of  each  particle  only  depends  on  the  overall 
force  acts  on  it  and  its  mass  [38].  In  the  GSA,  the  individuals  are 
a  collection  of  objects  also  the  solutions  of  the  problem  and  their 
corresponding  fitness  values  are  attributed  to  the  positions  and 
masses  of  those  objects.  The  GSA  expresses  the  gravitational  force 
as  follows: 


pk,t 

rje 


=  c‘ 


,  Mf  x  Mk 

k  x  — ^-j - ^x 

«/,+< 


(xf 


-x; 


k.t 


(29) 


where  the  gravitational  constant  decreases  during  the  optimization 
process  as: 


Ef  c  =  Y.  re  x  Eff  (37) 

eekbest,e=/=j 

where  kbest  is  a  set  of  the  particles  with  best  fitness  values.  The 
number  of  the  particles  in  kbest  is  a  function  of  time  which  starts 
with  Nswrm  and  decreases  linearly  to  1. 

Analyzing  the  GSA  algorithm,  the  following  points  can  be 
deduced: 

1.  In  GSA  algorithm,  it  is  expected  that  particles  are  attracted  by 
the  heaviest  (i.e.  most  optimal)  one  because  the  heavier  masses 
exert  more  powerful  gravitational  force  according  to  (29). 

2.  As  discussed  before,  according  to  the  Newton’s  law,  the  gravi¬ 
tational  force  between  two  particles  j  and  e  is  inversely 
proportional  to  the  (R/e)2  while  in  (29),  Rjeis  used  instead  of 
(Rje)2.  It  is  because  after  several  experiments  it  is  found  that 
this  displacement  provides  better  solutions  [29]. 

3.  Considering  (34),  the  inertia  mass  tends  to  make  the  mass 
movement  slow,  so  the  motion  of  the  heavier  masses  corre¬ 
sponded  to  the  better  solutions  is  more  slow  than  the  lighter 
ones.  Thus,  the  GSA  algorithm  searches  the  space  around  the 
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optimal  solutions  more  carefully  which  enhances  the  exploi¬ 
tation  and  local  search  capability  of  the  algorithm. 

4.  The  gravitational  constant  is  employed  to  make  balance 
between  exploration  and  exploitation.  It  has  a  big  value  at  the 
beginning  to  improve  the  exploration  power  of  the  algorithm 
and  help  to  avoid  being  trapped  in  local  optima.  Also,  it 
decreases  during  the  optimization  process  to  search  the  space 
with  the  higher  probability  of  the  optimal  solution  occurrence 
accurately. 


*“/„(x) 


i, 


/rx  -f„(x) 

fmax  _  fmin  5 
Jn  Jn 


/n(X)</nmin 

fn(X)>frX 

/„min  </„(X)  </nmax 


(41) 


The  DM’s  preference  is  defined  by  weight  factors  which  are 
allocated  to  the  objectives  then  the  normalized  membership  value 
for  each  solution  is  calculated  as: 


5.2.  Mutation 

This  paper  proposes  a  mutation  technique  to  improve  the 
convergence  characteristics  of  the  original  GSA.  This  mutation 
technique  is  proposed  to  improve  the  diversity  of  the  solutions, 
alleviate  the  stagnation  and  avoid  being  trapped  in  local  optima. 
For  each  particle,  three  particles  are  selected  randomly  as 
J=n2=£n3  =^j,  and  a  trial  solution  is  created  as: 

Xtmil  =  Xnew,n,  +  r  X  (X^,„2  -  X^w,„3)  (38) 


Nfid) 


ES^lWn  XJi/„(X,-) 

Ef=  1  £n°=l  WnX/fcft) 


(42) 


The  number  of  the  Pareto-optimal  solutions  of  the  problem  may 
be  excessive  or  even  infinitive.  Large  size  of  the  Pareto-optimal 
solution  causes  computational  burden  increase.  Furthermore, 
there  is  always  memory  constraint.  This  paper  uses  the  clustering 
approach  to  decrease  the  repository  members  without  destroying 
its  characteristics  [39]. 


7.  Solution  methodology 


Using  the  following  scheme,  a  mutant  solution  is  achieved: 


f  v*,t 

yt  _  I  *  trial  jd 
xmutjd  ~  \  k.t 

{  AnewJ0 


if  (r-)  <  r2) 
else 


(39) 


where  6  =  1, 2, n  and  r\  and  r2  are  two  random  numbers  with 
normal  distribution  between  0  and  1.  Finally,  for  each  object, 
between  the  mutant  position  and  position  found  by  the  original 
GSA,  the  one  with  the  better  fitness  value  is  selected  to  participate 
in  the  next  generation. 


6.  Multi-objective  MGSA 

Multi-objective  optimization  algorithms  are  devised  to  solve 
optimization  problems  with  several  incommensurable  and  con¬ 
flicting  objectives  stimulatingly.  These  algorithms,  instead  of  one 
single  solution  as  optimal  solution,  provide  a  set  of  solutions  named 
Pareto-optimal  solutions.  In  absence  of  the  any  preferable  infor¬ 
mation,  each  solution  in  this  set  has  no  priority  over  the  others.  The 
Pareto-optimal  solutions  are  selected  from  the  non-dominated 
solutions.  To  know  the  non-domination  concept,  suppose  that  Xi 
and  X2  are  two  optimal  solutions  for  a  given  multi-objective 
problem,  Xi  dominates  X2  if  and  only  if  the  following  two  condi¬ 
tions  are  satisfied: 

V/e{l,2,  ...,n0},  /KX^/KXa) 

3 qe {1,2,..., no},  /,(Xt)  </,(X2)  ^ 


7.1.  The  procedure  for  implementing  the  probabilistic  multi¬ 
objective  MGSA  can  be  summarized  in  the  following  steps 

Step  1 :  Randomly  generate  the  initial  positions  of  the  particles 
in  the  feasible  range  (6). 

Step  3:  Implement  the  2m  +  1  scheme.  Calculate  the  first  and 
the  second  moments  of  the  total  operation  cost  and  emissions. 
Step  4:  Initialize  the  repository  by  storing  the  non-dominated 
solutions  of  the  initial  population. 

Step  5:  Sort  the  particles  based  on  their  first  moment  values. 
Thereafter,  determine  worst k  and  best k. 

Step  6:  Update  the  gravitational  mass  of  each  particle  using 
(31)  and  (32). 

Step  7:  Update  the  overall  force  which  is  applied  on  each  particle 
using  (37). 

Step  8:  Calculate  the  acceleration  and  the  velocity  of  each 
particle  using  (34)  and  (35). 

Step  9:  Update  the  particles  position  using  (36). 

Step  10:  Implement  the  2m  +  1  scheme  as  Step  3. 

Step  11 :  Check  for  non-domination.  Update  the  repository. 

Step  12:  Apply  the  mutation  approach.  Find  the  new  solution  for 
each  particle. 

Step  13:  Implement  the  2m  +  1  scheme  as  Step  3. 

Step  14:  Check  for  non-domination.  Update  the  repository. 

Step  15:  Go  to  Step  4  until  the  current  iteration  number  reaches 
the  pre-specified  maximum  iteration  number. 

8.  Simulation  results 


This  paper  employs  the  MGSA  and  develops  it  to  devise 
a  powerful  multi-objective  algorithm  to  solve  the  EEM  of  the  MGs 
with  two  conflicting  objectives  as  emission  and  cost  minimization. 
To  do  this,  several  metaheuristic  techniques  are  used  to  find 
uniformly  POF  including  extreme  points  of  the  trade  of  surface. 
Firstly,  an  external  archive,  named  Repository,  is  planned  to  store 
the  non-dominated  solutions.  Each  solution  found  by  the  algorithm 
is  compared  with  the  other  ones  also  the  repository  members  for 
non-domination  and  the  non-dominated  solutions  are  stored  in  the 
repository. 

In  order  to  determine  the  best  compromise  solution,  a  fuzzy 
approach  is  used  in  this  paper.  In  the  fuzzy  approach,  a  fuzzy 
membership  function  is  allocated  to  each  of  the  objective  functions 
as  follows: 


The  MG  considered  in  this  paper  consists  of  different  DG  units 
such  as  an  MT,  FCPP,  PV,  and  WT.  It  is  supposed  that  all  DG  units 
produce  active  power  at  unity  power  factor,  neither  requesting  nor 
producing  reactive  power.  Furthermore,  there  is  a  power  exchange 
link  between  the  mentioned  MG  and  the  utility  (LV  network)  in 
order  to  trade  energy  during  a  day  based  on  decisions  of  the  Micro 
Grid  Central  Controller  (MGCC). 

Table  1  offers  different  parameters  used  in  this  paper  [4,15,30]. 
Parameters  of  the  FC  stack  presented  in  Table  1  are  adopted  from 
Ref.  [4].  The  prices  of  DER  units  are  considered  as  presented  in 
Ref.  [30]  while  the  prices  related  to  FCPP  are  estimated  by  scaling 
data  given  in  Ref.  [15].  The  normalized  forecasted  output  power  of 
WT  and  PV  units  and  the  hourly  forecasted  market  prices  for 
a  typical  day  can  be  found  in  Refs.  [30,40].  The  hourly  electrical  and 
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Table  1 

Simulation  parameters. 

Parameter 

Value 

T 

333  K 

A 

64  cm2 

Im 

178  pm 

Pu2 

1  (atm) 

Po2 

0.2095  (atm) 

B 

0.016 

Rc 

0.0003 

-0.948 

£2 

0.00286  +  0.0002  ln(A)  +  ln(Ctt2) 

h 

7.6  x  10  5 

U 

-1.93  x  10"4 

J  max 

469  mA  cm-2 

Jmax 

30  A 

Min,  max  power  of  FCPP 

10,  250  (kW) 

Min,  max  power  of  MT 

10,  80  (kW) 

Min,  max  power  of  PV 

0,  45  (kW) 

Min,  max  power  of  WT 

0,  35  (kW) 

Min,  max  power  of  Utility 

-95,  95  (kW) 

Price  of  natural  gas  for  FCPP 

0.21  (€ct/kWh) 

Fuel  price  for  residential  loads 

0.29  (€ct/kWh) 

Hydrogen  selling  price 

9.44  (€ct/l<Wh) 

Hydrogen  storage  efficiency 

0.95(€ct/kg) 

Hydrogen  pumping  cost 

0.09  (€ct/l<Wh) 

The  fuel  cell  cooling  time  constant 

0.75 

Hot  start  up  cost 

0.26  (€ct) 

Cold  start  up  cost 

0.7  (€ct) 

Bid  of  MT 

0.457  (€ct/kWh) 

Bid  of  PV  unit 

2.584  (€ct/l<Wh) 

Bid  of  WT 

1.073  (€ct/l<Wh) 

Nswrm 

40 

Iterm  ax 

300 

thermal  load  profiles  are  adopted  from  Ref.  [15].  Furthermore,  if  the 
recovered  thermal  energy  from  the  reformer  is  not  enough  to 
satisfy  the  thermal  load  (cooling,  heating,  or  hot  water)  additional 
heat  has  to  be  provided  by  an  auxiliary  boiler  of  the  CHP  system. 
The  data  for  the  auxiliary  boiler  can  be  found  in  Ref.  [31  ].  Finally, 
Table  2  presents  emission  coefficients  of  different  DER  units  [30,41  ]. 

In  the  following  subsections,  firstly,  in  order  to  show  the  effects 
of  hydrogen  production  in  the  EEM  problem,  deterministic  analysis 
are  performed  under  two  scenarios.  In  the  first  scenario  Si,  no 
hydrogen  is  produced  while  in  the  second  one  S2,  hydrogen 
production  strategy  is  employed.  The  multi-objective  EEM  problem 
is  solved  for  both  mentioned  scenarios  and  the  results  in  terms  of 
the  both  cost  and  emission  reductions  are  analyzed  in  detail. 
Finally,  the  probabilistic  multi-objective  MGSA  is  employed  to  find 
the  POF  of  the  probabilistic  EEM  problem  with  the  hydrogen 
production  strategy. 

8.1.  Deterministic  analysis 

In  the  deterministic  EEM  of  the  MGs,  it  is  assumed  that  the 
output  powers  of  the  WT  and  the  PV  units  are  equal  to  their  fore¬ 
casted  values  and  the  remaining  part  of  the  load  demand  is  satisfied 
by  the  other  DG  units.  Moreover,  market  prices  in  each  interval, 
natural  gas  price  for  the  FCPP,  operation  temperature  of  the  FC, 


Table  2 

Emission  coefficient  of  the  DG  sources. 


ID 

Type 

C02  (kg/MWh) 

S02  (kg/MWh) 

NO*  (kg/MWh) 

1 

MT 

720 

0.0036 

0.1 

2 

FCPP 

502.4780 

0.0036 

0.5215 

3 

PV 

0 

0 

0 

4 

WT 

0 

0 

0 

6 

Utility 

921.0585 

3.5827 

2.2947 

Cost  (Ect) 

Fig.  3.  POF  of  the  deterministic  EEM  problem  for  both  scenarios. 

pressure  of  the  reactant  gases  of  the  FC,  fuel  price  for  the  local 
thermal  loads  which  are  not  satisfied  by  the  FCPP  and  hydrogen 
selling  price  are  considered  to  be  equal  to  their  expected  values. 

Deterministic  multi-objective  MGSA  is  used  to  find  the  POF  of  the 
EEM  problem  of  the  mentioned  MG  under  aforesaid  scenarios.  The 
obtained  POF  for  both  S 1  and  S2  along  with  their  best  compromise 
solutions  are  portrayed  in  Fig.  3.  It  should  be  noted  that  in  finding 
best  compromise  solutions,  it  is  assumed  that  the  DM’s  preference  is 
unbiased,  that  is  to  say  both  wi  and  w2  are  equal  to  0.5.  As  shown 
Fig.  3,  great  improvement  in  terms  of  the  both  cost  and  emission 
reduction  using  the  hydrogen  production  and  storage  strategy  is 
evident.  The  POF  of  the  S2  is  well-distributed  in  the  cost-emission 
space.  Referring  Fig.  3,  for  equal  cost,  S2gives  lower  emission  in 
comparison  with  Si,  at  the  same  time,  for  equal  emission  for  both  Si 
and  S2,  S2  provides  lower  amount  of  emissions.  As  can  be  seen  in 
Fig.  3,  the  superiority  of  the  S2  over  Si  is  more  manifest  in  cost 
minimization.  For  more  convenience,  Table  3  tabulates  the  extreme 
points  of  the  trade-off  surface  also  the  best  compromise  solutions  for 
both  Si  and  S2.  By  analyzing  Table  3,  hydrogen  utilization  manage¬ 
ment  can  save  $621.7  in  per  day.  It  is  worthwhile  noting  that  the 
thermal  and  electrical  load  profiles  are  given  as  spring/fall  loads  in 
Ref.  [15],  so  the  proposed  hydrogen  management  can  save  about 
$111,906  in  spring  and  fall  seasons  which  is  surely  noteworthy.  In 
addition,  when  emission  minimization  is  the  goal  of  the  MGCC, 
hydrogen  management  saves  100  kg  daily  and  so,  18,000  kg  for  both 
spring  and  fall  seasons.  When  the  DM’s  preference  is  unbiased  S2 
saves  $439.7  while  its  produced  emission  is  just  114.26  kg  more  than 
Si.  Moreover,  with  the  same  cost  ($1466.84)  for  both  Si  and  S2,  well 
hydrogen  production  management  can  reduce  60  kg  emission 
production.  For  further  discussion  about  how  hydrogen  utilization 
strategy  makes  remarkable  profit  for  MGs,  Figs.  4-9  are  offered. 
Fig.  4  portrays  thermal  load  of  the  MG,  recovered  thermal  power  of 
the  reformer  and  generated  power  of  the  DG  units  for  the  power 
management  corresponding  to  the  minimum  cost  using  Si.  Fig.  5 
shows  the  same  data  for  the  power  management  corresponding  to 
the  minimum  emission  (Si ).  Moreover,  Fig.  6  gives  the  hourly  data  of 


Table  3 

Comparison  of  operational  cost  and  emission  of  scenarios  1  and  2. 


Method 

Best  cost 

Best 

emission 

Compromise 

solution 

Scenario  1 

Cost  (€ct) 

1179.29 

2518.74 

1466.84 

Emission  (kg) 

2391.08 

1889.69 

2.16952 

Scenario  2 

Cost  (€ct) 

557.53 

2418.63 

1027.14 

Emission  (kg) 

3083.46 

1844.24 

2283.78 
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Fig.  4.  Minimum  cost  -  (a)  thermal  load  of  the  MG  and  recovered  thermal  power,  (b)  generated  power  of  DER  units. 


Scenario  1  3 


Fig.  5.  Minimum  emission  -  (a)  thermal  load  of  the  MG  and 


recovered  thermal  power,  (b)  generated  power  of  DER  units. 


thermal  load  power  and  thermal  generated  power  by  the  reformer 
(Fig.  6a),  hydrogen  equivalent  production  (Fig.  6b),  the  secondary 
hydrogen  stream  amount  (Fig.  6c)  and  the  stored  hydrogen  (Fig.  6d) 
over  a  24  h  period  for  the  power  and  hydrogen  management 


corresponding  to  minimum  cost  (52).  In  addition,  Fig.  7  demon¬ 
strates  the  same  data  for  the  management  with  the  minimum 
emission  production.  Figs.  8  and  9  show  the  generated  power  of  the 
DG  units  for  minimum  cost  and  emission,  respectively,  when  S2  is 


Time(h) 


b 


Time(h) 


Time(h)  Time(h) 


Fig.  6.  Minimum  cost  —  (a)  thermal  load  and  recovered,  (b)  Hydrogen  production,  (c)  Secondary  hydrogen  stream  amount,  (d)  Hydrogen  equivalent  storage. 
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Fig.  7.  Minimum  emission  -  (a)  thermal  load  and  recovered,  (b)  Hydrogen  production,  (c)  Secondary  hydrogen  stream  amount,  (d)  Hydrogen  equivalent  storage. 


employed.  The  idea  behind  the  hydrogen  management  strategy  is 
storing  hydrogen  in  high  thermal  load  intervals  and  using  it  along 
with  hydrogen  produced  by  the  reformer  to  generate  electrical 
power  in  low  thermal  demand  intervals.  Figs.  6  and  7  in  comparison 
with  Figs.  4  and  5  clearly  show  that  the  proposed  methodology 
successfully  follows  the  aforementioned  strategy.  As  can  be  seen  in 
the  same  figures,  the  MG  experiences  low  thermal  load  between  14 
pm  and  18  pm.  In  Si,  during  14  pm-18  pm,  in  order  to  satisfy  elec¬ 
trical  load,  the  reformer  is  forced  to  produce  hydrogen  for  the  FCPP. 
Producing  hydrogen  is  followed  by  producing  heat  which  is  more 
than  the  heating  requirements  of  the  MG,  so  the  surplus  of  the 
thermal  power  to  the  MG’s  thermal  requirement  will  be  wasted.  This 
fact  causes  lose  money  and  generating  more  emissions.  On  the  other 
hand,  by  handling  hydrogen  utilization,  the  MGCC  can  achieve  more 
profitable  management.  As  Figs.  6  and  7  demonstrate,  in  the  first 
hours  of  the  day  which  MG  experiences  low  electrical  load,  the 
reformer  is  forced  to  produce  more  hydrogen  than  is  required  for 
FCPP  and  store  it  in  the  hydrogen  tank.  The  stored  hydrogen  is  used 
in  the  low  thermal  load  intervals,  more  between  14  pm  and  18  pm,  to 
supply  electrical  load  demand.  By  this  policy,  in  low  thermal  load 
intervals  the  reformer  doesn’t  have  to  produce  hydrogen  near  to  its 


Time(h) 


Fig.  8.  Minimum  cost  -  generated  electrical  power  of  different  DER  units. 


rated  capacity  to  supply  electrical  load  which  causes  to  waste  some 
amounts  of  its  generated  heat.  Moreover,  Fig.  8  in  comparison  with 
Fig.  4  shows  that  hydrogen  management  can  provide  more  profit¬ 
able  strategy  for  participating  in  open  market.  The  MGCC  buys 
power  during  low-price  periods  of  the  market  and  uses  the  capacity 
of  the  reformer  to  produce  hydrogen  for  storage.  This  stored 
hydrogen  can  be  used  to  produce  electrical  power  and  be  sold  to  the 
market  in  high-price  intervals.  This  strategy  can  make  a  good  profit 
for  the  MG  especially  in  reducing  operational  cost. 

In  order  to  verify  impacts  of  the  hydrogen  management  and 
thermal  recovery  on  the  efficiency  of  the  FCPP  Fig.  10  is  given.  To 
clarify  the  point  two  cases  are  defined.  In  the  first  case,  neither 
hydrogen  production  management  nor  thermal  energy  recovery  is 
considered  while  in  the  second  one  both  of  them  are  considered.  As 
shown  Fig.  10,  hydrogen  production  management  and  thermal 
power  recovery  improve  the  overall  efficiency  of  the  FCPP  greatly. 

In  order  to  make  better  sense  of  what  is  happening  Table  4  is 
given.  The  detailed  information  and  power  balances  of  some  hours 
are  tabulated  in  this  table.  This  information  are  regarding  to  best 


Fig.  9.  Minimum  emission  -  generated  electrical  power  of  different  DER  units. 
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cost  considering  S2.  The  design  parameters  of  fuel  cell  are  consid¬ 
ered  A  and  iFC.  A  is  an  adjustable  parameter  which  shows  the 
membrane  humidity.  This  parameter  depends  on  the  preparation 
procedure  of  the  membrane  and  can  be  regulated  by  the  relative 
humidity  and  Stoichiometry  ratio  of  the  anode  feed  gas.  iFc  is  cell 
operation  current.  Once  A  and  iFc  are  explored  based  on  the 
formulations  given  in  Ref.  [4],  Pceii  can  be  calculated.  For  the  other 
DGs  active  power  is  considered  as  design  variables.  Analyzing  this 
table  shows  how  electrical  and  thermal  load  are  fulfilled.  Further¬ 
more,  Fig.  11  shows  the  electrical  load  and  generation  for  all 
members  of  the  POFs  portrayed  in  Fig.  3.  As  shown  this  figure 
electrical  load  demands  of  all  hours  are  fulfilled  satisfactorily. 


8.2.  Probabilistic  analysis 

In  the  deterministic  analysis,  it  is  assumed  that  all  IRVs  values 
are  equal  to  their  forecasted  values.  Flowever,  as  discussed  in  the 
pervious  sections  some  required  data  for  the  EEM  can’t  forecasted 
exactly.  For  example,  wind  and  solar  power  have  stochastic  nature. 
Moreover,  in  the  open  access  market,  there  are  always  errors  in 
forecasted  values  for  power  market  and  fuel  prices.  Likewise,  the 
load  demands  are  more  unpredictable  than  before.  As  a  result,  the 
solutions  of  the  deterministic  optimization  algorithms  discussed  in 
the  previous  subsection  are  not  reliable  in  the  new  circumstance. 

In  this  section,  the  2m  +  1  point  estimate  method  is  imple¬ 
mented  along  with  the  multi-objective  MGSA  to  find  the  POF  of  the 
probabilistic  EEM  problem.  The  data  considered  as  the  uncertain 
inputs  are  as:  electrical  and  thermal  load  demands,  available  output 
power  of  the  PV  and  WT  units,  market  prices  in  each  interval, 
natural  gas  price  for  the  FCPP,  operation  temperature  of  the  FC, 
pressure  of  the  reactant  gases  of  FC,  fuel  price  for  the  thermal  loads 
of  the  MG  which  are  not  satisfied  by  the  FCPP  and  hydrogen  selling 


Fig.  11.  Electrical  load  and  generation  for  all  simulations. 


price.  The  normal  distribution  is  assumed  for  all  IRVs  except  output 
power  of  WT  unit  which  is  modeled  by  the  Beta  PDF  [25,42]. 

Fig.  12  shows  the  POF  of  the  probabilistic  EEM  problem  found  by 
the  MGSA  linked  to  the  2m  + 1  point  estimate  method.  Referring  to 
Fig.  12,  the  proposed  approach  successfully  found  uniformly  POF. 
The  DM  based  on  his/her  preference  and  priorities  selects  one  of 
the  members  of  the  Pareto-optimal  solutions.  Figs.  13  and  14 
portray  the  cumulative  density  functions  (CDFs)  of  the  cost  and 
emission  for  three  POF  solutions  (marked  in  Fig.  12).  The  Gram- 
Charlier  is  employed  to  find  the  accurate  distribution  for  total 
operational  cost  and  emission.  Figs.  13  and  14  demonstrate  that  the 
proposed  probabilistic  approach  could  successfully  describe  the 
variations  of  the  cost  and  emission  which  are  the  consequence  of 
the  variations  of  the  IRVs.  It  also  helps  the  system  operators  to 
know  how  likely  uncertainties  affect  the  system  and  how  to  handle 
the  related  issues. 


Table  4 

Detailed  information  for  some  hours  regarding  to  best  cost  scenario  2. 


FCPP 

MT 

PV 

WT 

Market 

Boiler 

Recovered 
thermal  energy 

Thermal  load 

Electrical 

load 

iFc 

X 

P cell 

1st  hour 

1.21 

9.51 

9.94 

6.06 

0 

4.17 

91.84 

93.32 

6.68 

100 

112 

2nd  hour 

3.02 

18.96 

23.08 

7.34 

0 

4.17 

79.42 

29.83 

73.17 

103 

114 

21st  hour 

27.78 

10.17 

161.69 

79.88 

0 

3.03 

-94.81 

4.01 

155.99 

160 

168 

22nd  hour 

28.59 

13.81 

163.96 

76.68 

0 

3.03 

-85.77 

10.53 

159.47 

170 

160 

23rd  hour 

9.34 

9.74 

64.31 

7.60 

0 

2.14 

80.56 

129.53 

46.47 

176 

155 

24th  hour 

8.32 

10.43 

57.96 

6.00 

0 

1.44 

84.58 

133.78 

41.22 

175 

150 
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Solution  1  Solution  2  Solution  3 


Fig.  13.  CDF  of  the  expected  cost. 


Solution  1 


Solution  2 


Solution  3 


Fig.  14.  CDF  of  the  expected  emission. 


Finally,  to  verify  the  impact  of  the  considering  uncertainty  in  the 
decision  making,  the  value  of  the  stochastic  solution  (VSS)  [43]  is 
calculated.  The  VSS  represents  how  much  the  MGCC  would  be 
willing  to  spend  in  order  to  know  the  future  realization  of  the 
problem  and  cover  the  uncertainty  in  the  stochastic  processes 
under  consideration.  The  VSS  is  calculated  by  subtracting  the  ex¬ 
pected  cost  (output  of  the  deterministic  EEM  problem)  from  the 
output  of  the  stochastic  model.  The  VSS  is  calculated  18  (€ct)  in  this 
test. 

9.  Conclusion 

In  this  paper,  a  probabilistic  approach  for  modeling  and  solving 
the  EEM  problem  of  renewable  MGs  under  uncertain  environment 
was  proposed.  In  the  proposed  model,  the  hydrogen  production 
and  utilization  management,  the  thermal  energy  recovery  and 
trade  with  the  upstream  grid  were  employed  along  with  different 
DERs  to  satisfy  the  thermal  and  the  electrical  load  simultaneously. 
The  thermal  energy  produced  by  the  reformer  was  recovered  while 


the  thermal  energy  of  the  FC  stack  was  neglected.  This  paper 
emphasized  particularly  the  PEMFCs  and  organized  the  economic/ 
emission  model  based  on  an  electrochemical  model  of  them.  A 
multi-objective  MGSA  was  proposed  and  linked  to  the  2m  +  1  point 
estimate  method  to  find  uniformly  POF  of  the  EEM  problem.  In 
order  to  verify  the  impact  of  the  hydrogen  production  on  the  EEM 
problem  two  scenarios  were  defined.  One  of  them  considered  the 
hydrogen  management  while  the  other  one  neglected  it.  The 
simulation  results  showed  the  great  improvement  which  the 
hydrogen  production  and  utilization  management  could  offer. 
The  enhancements  in  both  cost  and  emission  minimization  were 
discussed.  In  addition,  the  results  showed  that  the  hydrogen 
management  and  thermal  recovery  could  provide  30-35%  increase 
in  the  efficiency  of  the  FCPP.  Moreover,  the  simulation  results 
demonstrated  that  the  optimization  algorithm  force  the  MG  to 
store  some  amount  of  produced  hydrogen  in  high  thermal  load 
intervals  and  use  it  to  generate  electricity  in  low  thermal  load 
intervals.  By  this  policy,  less  thermal  power  would  be  wasted,  so 
more  money  and  emission  could  be  saved.  Finally,  to  cover  different 


298 


T.  Niknam  et  al.  /  Journal  of  Power  Sources  229  (2013)  285-298 


uncertainties  in  the  EEM  problem,  the  IRVs  are  modeled  by  the  Beta 
and  normal  distributions  and  the  probabilistic  POF  of  the  EEM 
problem  was  found.  The  CDF  of  three  members  of  the  probabilistic 
POF  were  portrayed.  The  figures  showed  that  the  proposed  prob¬ 
abilistic  approach  could  successfully  describe  the  variations  of  the 
cost  and  emission  which  follow  the  variation  of  the  IRVs.  It  also 
helps  the  system  operators  to  know  how  likely  uncertainties  affect 
the  system  and  how  to  handle  the  related  issues. 
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